В данном задании вам необходимо реализовать алгоритм кластеризации Partition Around Medoids.
Два возможных варианта реализации:
Нужно написать функцию, которая принимает на вход несколько параметров и возвращает также несколько значений.
Параметры функции:
pdist());Возвращаемые значения:
По аналогии с классами в scikit-learn, нужно реализовать класс, наследуемый от Base Estimator.
Подробнее про реализацию своих моделей в scikit-learn: here.
Параметры:
pdist());Методы:
Атрибуты:
Note 1: Параметры max_iter и tol должны иметь дефолтные значения.
Note 2: Функции для вычисления расстояний самим реализовывать не нужно, используйте pdist().
Также необходимо написать документацию к функции/методу: описать формат входных данных (параметров) и возвращаемых значений, особенности работы функции и детали реализации алгоритма. В качестве образца можно взять документацию к функциям/методам, которые рассматривали на занятиях.
Наивная реализация алгоритма будет работать довольно медленно - это нормально. Будет плюсом (но не является обязательным), если вы попытаетесь оптимизировать ваш код. Можете указать все ваши решения для оптимизации в документации.
import numpy as np
from sklearn.preprocessing import normalize
import scipy.spatial.distance as dis
from sklearn.decomposition import PCA
def create_medoids(X, k):
n_samples, n_features = np.shape(X)
medoids = np.zeros((k, n_features))
for i in range(k):
medoid = X[np.random.choice(range(n_samples))]
medoids[i] = medoid
non_medoids = []
for sample in X:
if not sample in medoids:
non_medoids.append(sample)
return medoids, non_medoids
def assign(X, medoids, k): #присваивает сэмплы к ближайшим медоидам
clusters = [[] for x in range(k)]
for sample_i, sample in enumerate(X):
closest_i = None
closest_distance = float("inf")
for i, medoid in enumerate(medoids):
distance = dis.euclidean(sample, medoid)
if distance < closest_distance:
closest_i = i
closest_distance = distance
medoid_i = closest_i
clusters[medoid_i].append(sample_i)
return clusters
def calculate_dist(X, clusters, medoids):
dist = 0
for i, cluster in enumerate(clusters):
medoid = medoids[i]
for sample_i in cluster:
dist += dis.euclidean(X[sample_i], medoid)
return dist
def pam(X, k):
X = [list(x) for x in X.to_numpy()]
medoids, non_medoids = create_medoids(X, k)
clusters = assign(X, medoids, k)
dist = calculate_dist(X, clusters, medoids)
while True:
best_medoids = medoids
lowest_dist = dist
for medoid in medoids:
for sample in non_medoids:
new_medoids = medoids.copy()
new_medoids[medoids == medoid] = sample
new_clusters = assign(X, new_medoids, k)
new_dist = calculate_dist(X, new_clusters, new_medoids)
if new_dist < lowest_dist:
lowest_dist = new_dist
best_medoids = new_medoids
if lowest_dist < dist:
dist = lowest_dist
medoids = best_medoids
else:
break
final_clusters = assign(X, medoids, k)
labels = np.zeros(np.shape(X)[0])
for cluster_i in range(len(clusters)):
cluster = clusters[cluster_i]
for sample_i in cluster:
labels[sample_i] = cluster_i
return labels
Алгоритм имеет две фазы:
Фаза Build:
Фаза Swap:
Схема реализации последовательного алгоритма (Псевдокод алгоритма)
1 функция PAM(D, k, tmax=100):
2 # D - матрица расстояний, k - число кластеров, tmax - маскимальное число итераций
3 выполнить фазу BUILD, получить множество метоидов M и множество не-метоидов L
4 вычислить значение целевой функции F
5 для t = 0..tmax-1:
6 выполнить фазу SWAP, вычислить значение целевой функции F'
7 delta = F - F'
8 если delta > 0:
9 обновить множества M и L
10 F = F'
11 иначе:
12 выйти из цикла
13 вернуть М
daata = norm
labels = pam(daata, 6)
labels
В рамках данной лабораторной работы вам предлагается проанализировать набор данных по различным городам США. Каждый город характеризуется следующими признаками:
import pandas as pd
pd.set_option('display.max_rows', 20)
pd.set_option('display.max_colwidth', 100)
data_desc = pd.read_csv('Data_Description.txt', sep=':')
data_desc
Housing и Crime - наоборот.Population- статистический признак, не имеющий интерпретации как “лучше-хуже”.Place - уникальный идентификатор объекта (города), он не должен использоваться при кластеризации.Longitude и Latitude. Их также не следует использовать при кластеризации данных.data = pd.read_csv('Data.txt', sep=' ')
data
0. Выполните необходимую предобработку данных. Перед кластеризацией исключите из данных признаки Place, Long и Lat.
есть
norm = data
norm = norm.drop(['Place'], axis=1)
norm = norm.drop(['Long'], axis=1)
norm = norm.drop(['Lat'], axis=1)
norm
from sklearn import preprocessing
names = norm.columns
scaler = preprocessing.StandardScaler()
norm = scaler.fit_transform(norm)
norm = pd.DataFrame(norm, columns=names)
norm
Посмотрим на различные расстояния между объектами применительно к этим данным
import numpy as np
from scipy.spatial.distance import pdist
import pandas as pd
from matplotlib import pyplot as plt
d = {'euclidean': pdist(norm, 'euclidean'),
'cityblock': pdist(norm, 'cityblock'),
'minkowski_6': pdist(norm, 'minkowski', p=6),
'cosine': pdist(norm, 'cosine'),
'chebyshev': pdist(norm, 'chebyshev'),
'canberra': pdist(norm, 'canberra')}
D = pd.DataFrame(d)
D.shape
D.plot(kind="bar",
subplots=True,
layout=(3,2),
sharey=False,
figsize=(15, 15)
)
plt.show()
На глаз результаты для всех расстояний выглядят похожими кроме canberra и косинусного. Посчитаем корреляции для более точной оценки.
D_corr = D.corr().loc[['cityblock', 'euclidean', 'minkowski_6', 'chebyshev', 'canberra', 'cosine'],
['cityblock', 'euclidean', 'minkowski_6', 'chebyshev', 'canberra', 'cosine']
]
D_corr
import seaborn as sns
plt.figure(figsize=(12, 8))
sns.heatmap(D_corr, annot=True)
plt.show()
Действительно, canberra и косинус сильно отличается от остальных, другие метрики на наших данных оказались очень похожими.
from scipy.cluster.hierarchy import linkage, dendrogram
import matplotlib.pyplot as plt
Z = linkage(norm, method='complete')
plt.figure(figsize=(15, 5))
dendrogram(Z, color_threshold=0)
plt.show()
type(Z)
Z
plt.figure(figsize=(10, 5))
dendrogram(Z, color_threshold=12.5)
plt.axhline(y=12.5, c='k')
plt.show()
plt.figure(figsize=(20, 10))
dendrogram(Z, color_threshold=2.5)
plt.axhline(y=2.5, c='k')
plt.show()
plt.figure(figsize=(20, 10))
dendrogram(Z, color_threshold=7.8)
plt.axhline(y=7.8, c='k')
plt.show()
import numpy as np
def plot_elbow(Z, h=13, w=5):
plt.figure(figsize=(h, w))
plt.plot(np.array(range(1, np.shape(Z)[0]+1)),
Z[:,2][::-1],
marker='o')
plt.xlabel("Number of clusters")
plt.ylabel("Merge distance")
plt.show()
plot_elbow(Z)
Выбираем ту точку, в которой снижение графика замедляется сильнее всего - "сгиб локтя". На графике это примерно 6 кластеров.
После того, как определились с порогом и числом кластеров, нужно извлечь полученное разбиение. Для этого используем функцию fcluster(). Функция возвращает массив длины 𝑛 , содержащий для каждого объекта i номер его кластера.
from scipy.cluster.hierarchy import fcluster
# По числу кластеров
fcluster(Z, t = 6, criterion='maxclust')
# По порогу (высоте дерева)
fcluster(Z, t = 7.8, criterion='distance')
Теперь можно проверить, удалось ли нам получить разбиение
import seaborn as sns
data['cluster'] = fcluster(Z, t = 7.8, criterion='distance')
sns.pairplot(data, hue='cluster', plot_kws={'alpha':0.5}, vars=['Climate',
'HousingCost',
'HlthCare',
'Crime',
'Transp',
'Educ',
'Arts',
'Recreat',
'Econ',
'Pop'])
plt.show()
В рассмотренном примере мы использовали Евклидово расстояние между объектами. Как мы увидели в начале анализа, расстояния Минковского, Чебышева и Манхеттенских кварталов почти идеально коррелируют с Евклидовым, и разумно ожидать аналогичных результатов кластеризации при их использовании.
Canberra очень сильно отличается от остальных. Давайте посмотрим, что получится при его использовании.
Z = linkage(norm, method='complete', metric='canberra')
plt.figure(figsize=(10, 10))
dendrogram(Z, color_threshold=9.5)
plt.axhline(y=9.5, c='k')
plt.show()
Косинусное расстояние:
Z = linkage(norm, method='complete', metric='cosine')
plt.figure(figsize=(10, 10))
dendrogram(Z, color_threshold=1.6)
plt.axhline(y=1.6, c='k')
plt.show()
from sklearn.cluster import DBSCAN
cl_dbscan = DBSCAN(eps=0.185, min_samples=4, metric='cosine')
cl_dbscan.fit(norm)
np.shape(cl_dbscan.labels_)
Посмотрим на то, как распределились объекты на корневые, пограничные и выбросы. Видно, что почти все объекты оказались корневыми, т.е. плотность объектов довольно высокая при выбранных параметрах.
data['points'] = 'Reachable'
data.iloc[cl_dbscan.core_sample_indices_, 4] = 'Core'
data.loc[cl_dbscan.labels_ == -1, 'points'] = 'Outlier'
sns.pairplot(x_vars="Long",
y_vars="Lat",
hue='points',
data=data,
height=10)
plt.show()
Посмотрим теперь на результаты кластеризации. Мы получили 6 кластеров и много выбросов. В целом полученная кластеризация выглядит разумно. Скопление объектов слева разделяется на несколько кластеров из-за сильного разброса по глубине эпицентра.
pd.Series(cl_dbscan.labels_).value_counts()
data['cl_dbscan'] = cl_dbscan.labels_
sns.pairplot(x_vars="Long",
y_vars="Lat",
hue='cl_dbscan',
data=data,
height=10)
plt.show()
sns.pairplot(x_vars='cl_dbscan',
y_vars=['Climate', 'HousingCost', 'HlthCare', 'Transp', 'Educ', 'Arts', 'Recreat', 'Econ', 'Pop'],
data=data,
height=5,
plot_kws={'alpha':0.15})
plt.show()
Эвристика для выбора параметром
from sklearn.neighbors import NearestNeighbors
neigh = NearestNeighbors(n_neighbors=2)
nbrs = neigh.fit(norm)
distances, indices = nbrs.kneighbors(norm)
distances = np.sort(distances, axis=0)
distances = distances[:,1]
plt.plot(distances)
cl_dbscan = DBSCAN(eps=3, min_samples=227)
cl_dbscan.fit(norm)
data['points'] = 'Reachable'
data.iloc[cl_dbscan.core_sample_indices_, 4] = 'Core'
data.loc[cl_dbscan.labels_ == -1, 'points'] = 'Outlier'
sns.pairplot(x_vars="Long",
y_vars="Lat",
hue='points',
data=data,
height=10)
plt.show()
pd.Series(cl_dbscan.labels_).value_counts()
data['cl_dbscan'] = cl_dbscan.labels_
sns.pairplot(x_vars="Long",
y_vars="Lat",
hue='cl_dbscan',
data=data,
height=10)
plt.show()
У меня не получилось евклидовым, какое значение соседей я бы не выбрала(
from sklearn.cluster import KMeans
data
norm
Z = KMeans(n_clusters = 6, # число кластеров
init = 'random',
n_init = 1,
max_iter = 100,
random_state=15434,
)
Z.fit(norm)
# координты центров кластеров
Z.cluster_centers_
# метки кластеров для каждого объекта
Z.labels_
# Значение целевой функции: сумма квадратов расстояний от объекта до центра кластера,
# к которому он принаджлежит
Z.inertia_
Добавим к данным координаты центров калстеров, чтобы отобразить их на графике. Чтобы не изменять исходные данные, создадим копию данных и будем работать с ней.
import numpy as np
data_kmeans = norm.copy()
centers = np.zeros((6, 10))
centers[:, 0:10] = Z.cluster_centers_
centers = pd.DataFrame(centers, columns=data_kmeans.columns)
data_kmeans['cluster'] = Z.labels_.astype(str)
centers['cluster'] = np.array(['c0', 'c1', 'c2', 'c3', 'c4', 'c5']).astype(str)
centers
data_kmeans = data_kmeans.append(centers, ignore_index=True)
data_kmeans.head(20)
sns.pairplot(data_kmeans.sort_values('cluster'),
hue='cluster',
plot_kws={'alpha':0.5},
vars=['Climate', 'HousingCost', 'HlthCare', 'Crime', 'Transp',
'Educ', 'Arts', 'Recreat', 'Econ', 'Pop']
)
plt.show()
inertia = []
for k in range(1, norm.shape[0]+1):
Z = KMeans(n_clusters=k,
init = 'random',
n_init = 100,
max_iter = 1000).fit(norm)
inertia.append(Z.inertia_)
plt.plot(range(1, data.shape[0]+1), inertia, 'bo-', marker='s')
plt.xlabel('$k$')
plt.ylabel('Objective function value')
plt.axvline(x=23, c='r')
plt.show()
inertia = np.sqrt(inertia)
plt.plot(range(1, norm.shape[0]+1), inertia, 'bo-', marker='s')
plt.xlabel('$k$')
plt.ylabel('$J(C_k)$')
plt.axvline(x=6, c='r')
plt.show()
При таких данных достаточно сложно определить точное количнство кластеров, я остановилась на 6
data_pam['cluster_Pam'] = pam(norm, 6)
plt.scatter(norm.to_numpy()[:, 0], norm.to_numpy()[:, 1],
c = [int(x) for x in data_pam['cluster_Pam'].to_numpy()],
s = 50, cmap = 'viridis')
!pip install hdbscan
import hdbscan
cl_hdbscan = hdbscan.HDBSCAN(min_cluster_size=2, gen_min_span_tree=True)
cl_hdbscan.fit(norm)
pd.Series(cl_hdbscan.labels_).value_counts()
import seaborn as sns
import matplotlib.pyplot as plt
data['cl_hdbscan'] = cl_hdbscan.labels_
sns.pairplot(x_vars="Long",
y_vars="Lat",
hue='cl_hdbscan',
data=data,
height=10)
plt.show()
plt.figure(figsize=(12, 8))
cl_hdbscan.minimum_spanning_tree_.plot(edge_cmap='viridis',
edge_alpha=0.6,
node_size=1,
edge_linewidth=2)
plt.show()
С помощью hdbscan удалось разделить максимум на 5 кластеров, что было близко к моим ассуждениям. Подбирая другие параметры, распределение было намного хуже, это - лучший вариант.
В результате выполнения предыдущих пунктов вы должны получить 4 или больше разбиения объектов (по одному на каждый метод). Сравните их между собой, сделайте выводы о сходствах и различиях.
Сравнивая различные методы кластеризации, я нашла несколько сходств и различий.
Что касается сходств, кластеров получилось примерно одно и то же количество(5 штук в hdbscan, 6 штук - в остальных). По правде говоря, я "подгоняла" значение 6 для различных кластеризаций, так как посчитала его оптимально верным. Также для меня достаточно ярко почти во всех методах выделился "островочек" точек, находящийся слева чуть ниже середины (например, выделен зеленым цветом в dbscan'е). Да и вообще, кластеризация в целом везде выглядит в общем похоже.
С другой стороны, естсественно, мы найдем и несколько различий. Во-первых, очевидно, что мы не получили абсолютно идентичную реализацию во всех методах. Во-вторых, время выполнения тоже отличается. Далее см. таблицу ниже.
from sklearn import metrics
from sklearn import datasets
import pandas as pd
from sklearn.cluster import KMeans, AgglomerativeClustering, AffinityPropagation, SpectralClustering
data = datasets.load_digits()
X, y = data.data, data.target
algorithms = []
algorithms.append(KMeans(n_clusters=6, random_state=1))
algorithms.append(AffinityPropagation())
algorithms.append(SpectralClustering(n_clusters=6, random_state=1,
affinity='nearest_neighbors'))
algorithms.append(AgglomerativeClustering(n_clusters=6))
data = []
for algo in algorithms:
algo.fit(X)
data.append(({
'ARI': metrics.adjusted_rand_score(y, algo.labels_),
'AMI': metrics.adjusted_mutual_info_score(y, algo.labels_),
'Homogenity': metrics.homogeneity_score(y, algo.labels_),
'Completeness': metrics.completeness_score(y, algo.labels_),
'V-measure': metrics.v_measure_score(y, algo.labels_),
'Silhouette': metrics.silhouette_score(X, algo.labels_)}))
results = pd.DataFrame(data=data, columns=['ARI', 'AMI', 'Homogenity',
'Completeness', 'V-measure',
'Silhouette'],
index=['K-means', 'Affinity',
'Spectral', 'Agglomerative'])
results
Выберите одно разбиение, наиболее подходящее на ваш взгляд. Предложите интерпретацию полученным кластерам или покажите, что этого сделать нельзя.
copy = data
data_null = copy
data_null = data_null.loc[data_null['cluster'] != '1']
data_null = data_null.loc[data_null['cluster'] != '2']
data_null = data_null.loc[data_null['cluster'] != '3']
data_null = data_null.loc[data_null['cluster'] != '4']
data_null = data_null.loc[data_null['cluster'] != '5']
data_null.head(20)
data_null.tail(20)
Анализируя кусочек данного датасета, можно заметить, что особо данные не схожи, так как их очень много. Но сходства можно увидеть, например, в колонках Educ, Econ, Transp.
copy = data
data_one = copy
data_one = data_one.loc[data_one['cluster'] != '0']
data_one = data_one.loc[data_one['cluster'] != '2']
data_one = data_one.loc[data_one['cluster'] != '3']
data_one = data_one.loc[data_one['cluster'] != '4']
data_one = data_one.loc[data_one['cluster'] != '5']
data_one.head(20)
data_one.tail(20)
Похожие значения в признаках Crime, Educ, Econ.
copy = data
data_two = copy
data_two = data_two.loc[data_two['cluster'] != '0']
data_two = data_two.loc[data_two['cluster'] != '1']
data_two = data_two.loc[data_two['cluster'] != '3']
data_two = data_two.loc[data_two['cluster'] != '4']
data_two = data_two.loc[data_two['cluster'] != '5']
data_two.head(20)
data_two.tail(20)
Все те же Educ и Econ, также вижу схожие показатели HousingCost.
copy = data
data_three = copy
data_three = data_three.loc[data_three['cluster'] != '0']
data_three = data_three.loc[data_three['cluster'] != '2']
data_three = data_three.loc[data_three['cluster'] != '1']
data_three = data_three.loc[data_three['cluster'] != '4']
data_three = data_three.loc[data_three['cluster'] != '5']
data_three.head(20)
data_three.tail(20)
Educ имеет примерно одно и то же значние, разброс невелик. В принципе, показатель Econ не сильно отличается в каждой из строчек. Да, разброс есть, но не сильно большой.
copy = data
data_four = copy
data_four = data_four.loc[data_four['cluster'] != '0']
data_four = data_four.loc[data_four['cluster'] != '2']
data_four = data_four.loc[data_four['cluster'] != '1']
data_four = data_four.loc[data_four['cluster'] != '3']
data_four = data_four.loc[data_four['cluster'] != '5']
data_four.head(20)
data_four.tail(20)
Educ имеет примерно одно и то же значние, разброс невелик. В принципе, показатель Econ не сильно отличается в каждой из строчек. Да, разброс есть, но не сильно большой.
copy = data
data_five = copy
data_five = data_five.loc[data_five['cluster'] != '0']
data_five = data_five.loc[data_five['cluster'] != '2']
data_five = data_five.loc[data_five['cluster'] != '1']
data_five = data_five.loc[data_five['cluster'] != '3']
data_five = data_five.loc[data_five['cluster'] != '4']
data_five
Схожие показатели HousingCost, Crime, Transp, Educ, Recreat, Econ. Неудивительно, что здесь нашлось больше схожих признаков, ведь у нас всего 2 строчки, вероятность совпадения выше.
У всех кластеров сошлись примерно показатель Educ и Econ, но также это из-за того, что вообще у всех кластеров они схожи, у значений этого признака разброс небольшой. Что же касается разделения на кластеры, они все имеют по одному схожему признаку: где-то Crime, где-то HousingCost и тд. Не знаю точно, но это мое предположение, почему так поделились данные. Про географическое расположение см ниже.
Оцените, как полученные кластеры распределены географически. (Бонусное) Провизуализируйте распределение на карте США. Оцените, как полученные кластеры распределены по штатам. Можно ли выделить какую-то зависимость (территориальную или для штатов)?
!pip install folium
import folium
data['cluster'] = data_kmeans['cluster']
data
m = folium.Map(
location=[45.372, -121.6972],
zoom_start=4,
tiles='Stamen Terrain'
)
for i, row in data.iterrows():
if row['cluster'] == '3':
color='red'
elif row['cluster'] == '1':
color='green'
elif row['cluster'] == '2':
color='white'
elif row['cluster'] == '0':
color='pink'
elif row['cluster'] == '4':
color='blue'
else: color='orange'
folium.Marker(location=[row['Lat'], row['Long']], popup=i, icon=folium.Icon(color=color)).add_to(m)
m
Как можно заметить, большинство данных было собрано с восточной части Соединенных Штатов (преобладает розовый, белый и красный кластеры). Западная половина немного пустует, но ближе к "берегам" там тоже сосредоточены кластеры, в особенности синего цвета.
Красный кластер в оснвном расположился в 4 четверти (юго-восток) и все точки располагаются в нижней части.
В целом зависимости четок не обнаружила, kmeans распределил по своему алгоритму. Мое предположение см. выше.